#GOALS
#1.Describe how to built UPSS.
#2.Exploring UPSS effects on canopy microclimate, including temperature and light intensity, relative to the open and shrub.
#3.Understanding how different light permeabilities and shelter shapes influence the above parameters.
library (tidyverse)
## Warning: package 'tidyverse' was built under R version 3.6.2
## -- Attaching packages ----------------------- tidyverse 1.3.0 --
## v ggplot2 3.2.1     v purrr   0.3.3
## v tibble  2.1.3     v dplyr   0.8.3
## v tidyr   1.0.0     v stringr 1.4.0
## v readr   1.3.1     v forcats 0.4.0
## Warning: package 'forcats' was built under R version 3.6.2
## -- Conflicts -------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(ggplot2)
#import datasets
panoche.climate.2019<-read.csv("C:/Users/Nargol Ghazian/Desktop/Animal-Behaviour-and-Climate-project/CH3/panoche_climate_hourly 2019.csv")

mean.daily.panoche <- read.csv("C:/Users/Nargol Ghazian/Desktop/Animal-Behaviour-and-Climate-project/CH3/weather station climate.csv")

str(panoche.climate.2019)
## 'data.frame':    8760 obs. of  28 variables:
##  $ station.id  : int  56 56 56 56 56 56 56 56 56 56 ...
##  $ site        : Factor w/ 2 levels "","panoche": 2 2 2 2 2 2 2 2 2 2 ...
##  $ station.name: Factor w/ 2 levels "","Los Banos": 2 2 2 2 2 2 2 2 2 2 ...
##  $ CIMIS.Region: Factor w/ 2 levels "","San Joaquin Valley": 2 2 2 2 2 2 2 2 2 2 ...
##  $ hour        : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ date        : Factor w/ 27 levels "","5/18/2019",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ month       : int  5 5 5 5 5 5 5 5 5 5 ...
##  $ day         : int  18 18 18 18 18 18 18 18 18 18 ...
##  $ year        : int  2019 2019 2019 2019 2019 2019 2019 2019 2019 2019 ...
##  $ precip      : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ radiation   : int  NA NA NA NA NA 57 176 432 597 723 ...
##  $ air.temp    : num  8.7 9.2 8.8 7.6 6.6 7 10.3 13.7 16 17.5 ...
##  $ wind.speed  : num  1.4 1.3 0.9 1.2 1.2 0.6 1 0.9 1.4 3 ...
##  $ soil.temp   : num  18.1 18 18 17.9 17.8 17.7 17.6 17.5 17.4 17.3 ...
##  $ air.temp.F  : num  47.7 48.6 47.8 45.7 43.9 ...
##  $ X           : logi  NA NA NA NA NA NA ...
##  $ X.1         : logi  NA NA NA NA NA NA ...
##  $ X.2         : logi  NA NA NA NA NA NA ...
##  $ X.3         : logi  NA NA NA NA NA NA ...
##  $ X.4         : logi  NA NA NA NA NA NA ...
##  $ X.5         : logi  NA NA NA NA NA NA ...
##  $ X.6         : logi  NA NA NA NA NA NA ...
##  $ X.7         : logi  NA NA NA NA NA NA ...
##  $ X.8         : logi  NA NA NA NA NA NA ...
##  $ X.9         : logi  NA NA NA NA NA NA ...
##  $ X.10        : logi  NA NA NA NA NA NA ...
##  $ X.11        : logi  NA NA NA NA NA NA ...
##  $ X.12        : Factor w/ 3 levels " ","R","Y": 1 1 1 1 1 1 1 1 1 1 ...
na.omit (panoche.climate.2019)#exclude missing values
##  [1] station.id   site         station.name CIMIS.Region hour        
##  [6] date         month        day          year         precip      
## [11] radiation    air.temp     wind.speed   soil.temp    air.temp.F  
## [16] X            X.1          X.2          X.3          X.4         
## [21] X.5          X.6          X.7          X.8          X.9         
## [26] X.10         X.11         X.12        
## <0 rows> (or 0-length row.names)
head(is.na(panoche.climate.2019))#check for missing values
##      station.id  site station.name CIMIS.Region  hour  date month   day
## [1,]      FALSE FALSE        FALSE        FALSE FALSE FALSE FALSE FALSE
## [2,]      FALSE FALSE        FALSE        FALSE FALSE FALSE FALSE FALSE
## [3,]      FALSE FALSE        FALSE        FALSE FALSE FALSE FALSE FALSE
## [4,]      FALSE FALSE        FALSE        FALSE FALSE FALSE FALSE FALSE
## [5,]      FALSE FALSE        FALSE        FALSE FALSE FALSE FALSE FALSE
## [6,]      FALSE FALSE        FALSE        FALSE FALSE FALSE FALSE FALSE
##       year precip radiation air.temp wind.speed soil.temp air.temp.F    X
## [1,] FALSE  FALSE      TRUE    FALSE      FALSE     FALSE      FALSE TRUE
## [2,] FALSE  FALSE      TRUE    FALSE      FALSE     FALSE      FALSE TRUE
## [3,] FALSE  FALSE      TRUE    FALSE      FALSE     FALSE      FALSE TRUE
## [4,] FALSE  FALSE      TRUE    FALSE      FALSE     FALSE      FALSE TRUE
## [5,] FALSE  FALSE      TRUE    FALSE      FALSE     FALSE      FALSE TRUE
## [6,] FALSE  FALSE     FALSE    FALSE      FALSE     FALSE      FALSE TRUE
##       X.1  X.2  X.3  X.4  X.5  X.6  X.7  X.8  X.9 X.10 X.11  X.12
## [1,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE
## [2,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE
## [3,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE
## [4,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE
## [5,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE
## [6,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE
head(na.omit (mean.daily.panoche))
##   station.name    site  sensor      lat    long      date year month day
## 1    los banos panoche station -120.867 37.0563 5/18/2019 2019     5  18
## 2    los banos panoche station -120.867 37.0563 5/19/2019 2019     5  19
## 3    los banos panoche station -120.867 37.0563 5/20/2019 2019     5  20
## 4    los banos panoche station -120.867 37.0563 5/21/2019 2019     5  21
## 5    los banos panoche station -120.867 37.0563 5/22/2019 2019     5  22
## 6    los banos panoche station -120.867 37.0563 5/23/2019 2019     5  23
##   temp.max temp.min temp.mean
## 1     83.1     55.2     69.15
## 2     83.3     55.4     69.35
## 3     83.5     55.4     69.45
## 4     83.7     55.6     69.65
## 5     83.9     55.8     69.85
## 6     84.2     55.9     70.05
#DATA VIZ MACRO (WEATHER STATION) 
##macro-climate plots
ggplot (panoche.climate.2019, aes((hour), air.temp)) + geom_smooth() + xlab("Hour (0-24)") + ylab ("Temperature (°C)")+  theme_classic()+ theme(axis.text=element_text(size=12))+ theme(axis.text=element_text(size=12))#how air temperature changed during 24h period in 2019 during study period
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## Warning: Removed 8136 rows containing non-finite values (stat_smooth).

ggplot (panoche.climate.2019, aes((day), air.temp)) + geom_smooth() + xlab("Day") + ylab ("Temperature (°C)")+  theme_classic()+ theme(axis.text=element_text(size=12))+ theme(axis.text=element_text(size=12))#air temperature over the days of the study period
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## Warning: Removed 8136 rows containing non-finite values (stat_smooth).

ggplot (panoche.climate.2019, aes((hour), soil.temp)) + geom_smooth() + xlab("Hour (0-24)") + ylab ("Temperature (°C)")+  theme_classic()+ theme(axis.text=element_text(size=12))+ theme(axis.text=element_text(size=12))#how soil temperature changed during 24h period in 2019 during study period
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## Warning: Removed 8136 rows containing non-finite values (stat_smooth).

ggplot (panoche.climate.2019, aes((day), soil.temp)) + geom_smooth() + xlab("Day") + ylab ("Temperature (°C)")+  theme_classic()+ theme(axis.text=element_text(size=12))+ theme(axis.text=element_text(size=12))#how soil temperature changed during 24h period in 2019 during study period
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## Warning: Removed 8136 rows containing non-finite values (stat_smooth).

ggplot (panoche.climate.2019, aes((hour), radiation)) + geom_smooth() + xlab("Hour (0-24)") + ylab ("Solar Radiation (W/m²)")+  theme_classic()+ theme(axis.text=element_text(size=12))+ theme(axis.text=element_text(size=12))#how solar radiation changed during 24h period in 2019 during study period
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## Warning: Removed 8347 rows containing non-finite values (stat_smooth).

ggplot (panoche.climate.2019, aes((day), radiation)) + geom_smooth() + xlab("Day") + ylab ("Solar Radiation (W/m²)")+  theme_classic()+ theme(axis.text=element_text(size=12))+ theme(axis.text=element_text(size=12))#sunlight experinced over the the study period
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## Warning: Removed 8347 rows containing non-finite values (stat_smooth).

library (ggbeeswarm)
## Warning: package 'ggbeeswarm' was built under R version 3.6.3
ggplot (panoche.climate.2019, aes(as.factor(date), air.temp)) + geom_boxplot() + xlab("Date") + ylab ("Temperature (°C)")+ geom_quasirandom(alpha=0.05)+ theme_classic()+ theme(axis.text=element_text(size=12))+ theme(axis.text=element_text(size=12))+ theme(axis.text.x = element_text(angle = 90))+ geom_smooth(se=FALSE, color="black", aes(group=1))#daily air temperature averages
## Warning: Removed 8136 rows containing non-finite values (stat_boxplot).
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## Warning: Removed 8136 rows containing non-finite values (stat_smooth).
## Warning: Removed 8136 rows containing missing values
## (position_quasirandom).

#DATA VIZ MICRO
##micro-climate plots
##let's try some plots for shelter, shrub, and open

shelter.shrub.open <- read.csv("C:/Users/Nargol Ghazian/Desktop/Animal-Behaviour-and-Climate-project/CH3/shrub_contrast_final.csv")#import dataset

ggplot(shelter.shrub.open, aes((microsite), temp, fill=microsite)) + geom_boxplot() + xlab("Microsite") + ylab ("Temperature (°F)")+ theme_classic()+ theme(axis.text=element_text(size=12))+theme(axis.text.x = element_text(angle = 90))+theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())+ labs(fill = "Microsite")+ stat_summary(fun.y=mean, colour="black", geom="point", shape=18, size=3,show_guide = FALSE)#Boxplot temperaure by microsite
## Warning: `show_guide` has been deprecated. Please use `show.legend`
## instead.

ggplot(shelter.shrub.open, aes((microsite), intensity, fill=microsite)) + geom_boxplot() + xlab("Microsite") + ylab ("Solar Radiation (lum/ft²)")+ theme_classic()+ theme(axis.text=element_text(size=12))+theme(axis.text.x = element_text(angle = 90))+theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())+ labs(fill = "Microsite")+ stat_summary(fun.y=mean, colour="black", geom="point", shape=18, size=3,show_guide = FALSE)#Boxplot light intensity by microsite
## Warning: `show_guide` has been deprecated. Please use `show.legend`
## instead.
## Warning: Removed 9330 rows containing non-finite values (stat_boxplot).
## Warning: Removed 9330 rows containing non-finite values (stat_summary).

ggplot(shelter.shrub.open, aes((day), temp, color=microsite)) + geom_smooth()+ xlab("Day") + ylab ("Temperature (°F)")+ theme_classic()+ theme(axis.text=element_text(size=12))+stat_summary(fun.y=max, geom="point", size=2, aes(shape = microsite))+ theme(legend.position = "none") + facet_grid("microsite")
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

ggplot(shelter.shrub.open, aes((day), intensity, color=microsite)) + xlab("Day") + ylab ("Solar Radiation (lum/ft²)")+ theme_classic()+ theme(axis.text=element_text(size=12))+geom_smooth()+ stat_summary(fun.y=max, geom="point", size=2, aes(shape = microsite))+ theme(legend.position = "none") + facet_grid("microsite")
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## Warning: Removed 9330 rows containing non-finite values (stat_smooth).

## Warning: Removed 9330 rows containing non-finite values (stat_summary).

library(ggpubr)
## Loading required package: magrittr
## 
## Attaching package: 'magrittr'
## The following object is masked from 'package:purrr':
## 
##     set_names
## The following object is masked from 'package:tidyr':
## 
##     extract
ggline(shelter.shrub.open, x = "day", y = "temp", color = "microsite",
 add = "mean_ci", shape = "microsite", xlab = "Day", ylab = "Temperature (°F)", legend.title= "Microsite", legend="right")#dotted line graph of mean daily temperatures

ggline(shelter.shrub.open, x = "day", y = "intensity", color = "microsite",
 add = "mean_ci", shape = "microsite", legend.title= "Microsite", xlab = "Day", ylab = "Solar Radiation (lum/ft²)", legend="right")#dotted line graph of mean daily solar radiation 
## Warning: Removed 9330 rows containing non-finite values (stat_summary).

#visualizing mean, median, and mode with histograms
ggplot(shelter.shrub.open, aes(temp, fill = microsite)) +
         geom_histogram(binwidth = 5) +
  scale_fill_brewer(palette = "Set1")+   labs(fill = "", x = "Temperature (°F)", y = "Frequency")+theme_classic()+ theme(axis.text=element_text(size=12))+labs(fill = "Microsite")

library(ggpubr)
ggplot(shelter.shrub.open, aes(temp, fill = microsite)) +
         geom_histogram(binwidth = 5) +
  scale_fill_brewer(palette = "Set1")+   labs(fill = "", x = "Temperature (°F)", y = "Frequency")+theme_classic()+ theme(axis.text=element_text(size=12))+ stat_central_tendency(aes(color = microsite), type = "mean", color="green", linetype = 2)+ stat_central_tendency(type = "median", color = "blue", linetype = 2)+ facet_grid(~microsite)+theme(legend.position = "none")

ggplot(shelter.shrub.open, aes(temp, fill = microsite)) +
         geom_density() +
  scale_fill_brewer(palette = "Set1")+  labs(fill = "", x = "Temperature (°F)", y = "Frequency")+theme_classic()+ theme(axis.text=element_text(size=12))+labs(fill = "Microsite")

ggplot(shelter.shrub.open, aes(temp, fill = microsite)) +
         geom_density() +
  scale_fill_brewer(palette = "Set1")+   labs(fill = "", x = "Temperature (°F)", y = "Frequency")+theme_classic()+ theme(axis.text=element_text(size=12))+ stat_central_tendency(aes(color = microsite), type = "mean", color="green", linetype = 2)+ stat_central_tendency(type = "median", color = "blue", linetype = 2)+ facet_grid(~microsite)+theme(legend.position = "none")

#EDA
#explore normality 
#can't do Shapiro-Wilk dataset too large 
library(ggpubr)

ggqqplot(shelter.shrub.open$temp)#Data are positively(right) skewed, Gaussian  

hist(shelter.shrub.open$temp) 

ggqqplot(shelter.shrub.open$intensity)#Data follow a poisson distribution
## Warning: Removed 9330 rows containing non-finite values (stat_qq).
## Warning: Removed 9330 rows containing non-finite values (stat_qq_line).

## Warning: Removed 9330 rows containing non-finite values (stat_qq_line).

hist(shelter.shrub.open$intensity)

#Neither intensity nor temperature follow a normal distribution
#explore light and temperature relationship
macro.micro.contrast <- read.csv("C:/Users/Nargol Ghazian/Desktop/Animal-Behaviour-and-Climate-project/CH3/micro-macro-contrast.csv")

cor.test(macro.micro.contrast$temp, macro.micro.contrast$intensity, method = "kendall")#light and temperature are correlated p<0.0001, tau=0.2813
## 
##  Kendall's rank correlation tau
## 
## data:  macro.micro.contrast$temp and macro.micro.contrast$intensity
## z = 45.612, p-value < 2.2e-16
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
##       tau 
## 0.2659057
#scatterplot for the relationship
ggplot(macro.micro.contrast, aes(temp, intensity))+ geom_point()+ ylab("Light Intensity (lum/ft²)")+ xlab("Temperature (°F)")+geom_smooth(method='glm', se=FALSE)+theme_classic()#Perhaps too dense too look at
## Warning: Removed 9330 rows containing non-finite values (stat_smooth).
## Warning: Removed 9330 rows containing missing values (geom_point).

ggplot(macro.micro.contrast, aes(temp, intensity))+geom_point()+ ylab("Light Intensity (lum/ft²)")+ xlab("Temperature (°F)")+geom_smooth(method='glm', se=FALSE)+  facet_grid(~microsite)+ theme_classic()#faceted by microsite
## Warning: Removed 9330 rows containing non-finite values (stat_smooth).

## Warning: Removed 9330 rows containing missing values (geom_point).

ggplot(macro.micro.contrast, aes(temp, intensity))+geom_point()+ ylab("Light Intensity (lum/ft²)")+ xlab("Temperature (°F)")+geom_smooth(method='glm', se=FALSE)+  facet_grid(~cover.type)+ theme_classic()#faceted by cover type
## Warning: Removed 9330 rows containing non-finite values (stat_smooth).

## Warning: Removed 9330 rows containing missing values (geom_point).

#MODELS
##model hypotheses
##pair-wise comparison
library(sjPlot)
## Warning: package 'sjPlot' was built under R version 3.6.3
## Install package "strengejacke" from GitHub (`devtools::install_github("strengejacke/strengejacke")`) to load all sj-packages at once!
lm.temp <- glm(temp~as.factor(microsite), data = shelter.shrub.open, family="gaussian")
summary (lm.temp)
## 
## Call:
## glm(formula = temp ~ as.factor(microsite), family = "gaussian", 
##     data = shelter.shrub.open)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -39.543  -16.496   -5.442   13.796  110.218  
## 
## Coefficients:
##                                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                        73.6902     0.2214 332.819  < 2e-16 ***
## as.factor(microsite)shrub.ambient   0.2471     0.4179   0.591   0.5544    
## as.factor(microsite)shrub.soil      3.3369     0.4760   7.011 2.44e-12 ***
## as.factor(microsite)square         -1.0049     0.4414  -2.277   0.0228 *  
## as.factor(microsite)triangle       -3.1949     0.5171  -6.179 6.56e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 483.3702)
## 
##     Null deviance: 10668889  on 21959  degrees of freedom
## Residual deviance: 10612393  on 21955  degrees of freedom
## AIC: 198057
## 
## Number of Fisher Scoring iterations: 2
tab_model(lm.temp)
  temp
Predictors Estimates CI p
(Intercept) 73.69 73.26 – 74.12 <0.001
microsite [shrub.ambient] 0.25 -0.57 – 1.07 0.554
microsite [shrub.soil] 3.34 2.40 – 4.27 <0.001
microsite [square] -1.00 -1.87 – -0.14 0.023
microsite [triangle] -3.19 -4.21 – -2.18 <0.001
Observations 21960
R2 Nagelkerke 0.924 
library (emmeans)
## Welcome to emmeans.
## NOTE -- Important change from versions <= 1.41:
##     Indicator predictors are now treated as 2-level factors by default.
##     To revert to old behavior, use emm_options(cov.keep = character(0))
anova(lm.temp, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: temp
## 
## Terms added sequentially (first to last)
## 
## 
##                      Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                                 21959   10668889              
## as.factor(microsite)  4    56496     21955   10612393 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans1 <- emmeans(lm.temp, pairwise~microsite)

lm.intensity <- glm(intensity~as.factor(microsite), data=shelter.shrub.open, family = "quasipoisson")
summary (lm.intensity)
## 
## Call:
## glm(formula = intensity ~ as.factor(microsite), family = "quasipoisson", 
##     data = shelter.shrub.open)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -81.516  -56.469  -33.933    2.815  293.607  
## 
## Coefficients:
##                                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                        8.11118    0.01800  450.56   <2e-16 ***
## as.factor(microsite)shrub.ambient -0.71605    0.04520  -15.84   <2e-16 ***
## as.factor(microsite)shrub.soil    -4.58886    0.41459  -11.07   <2e-16 ***
## as.factor(microsite)square        -0.68732    0.04727  -14.54   <2e-16 ***
## as.factor(microsite)triangle      -0.58210    0.05431  -10.72   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 6402.7)
## 
##     Null deviance: 62739687  on 12629  degrees of freedom
## Residual deviance: 54663384  on 12625  degrees of freedom
##   (9330 observations deleted due to missingness)
## AIC: NA
## 
## Number of Fisher Scoring iterations: 6
tab_model(lm.intensity)
  intensity
Predictors Incidence Rate Ratios CI p
(Intercept) 3331.50 3215.33 – 3450.44 <0.001
microsite [shrub.ambient] 0.49 0.45 – 0.53 <0.001
microsite [shrub.soil] 0.01 0.00 – 0.02 <0.001
microsite [square] 0.50 0.46 – 0.55 <0.001
microsite [triangle] 0.56 0.50 – 0.62 <0.001
Observations 12630
R2 Nagelkerke 1.000 
anova(lm.intensity, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: quasipoisson, link: log
## 
## Response: intensity
## 
## Terms added sequentially (first to last)
## 
## 
##                      Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                                 12629   62739687              
## as.factor(microsite)  4  8076303     12625   54663384 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(lm.temp, pairwise~microsite)
## $emmeans
##  microsite     emmean    SE  df asymp.LCL asymp.UCL
##  open            73.7 0.221 Inf      73.3      74.1
##  shrub.ambient   73.9 0.354 Inf      73.2      74.6
##  shrub.soil      77.0 0.421 Inf      76.2      77.9
##  square          72.7 0.382 Inf      71.9      73.4
##  triangle        70.5 0.467 Inf      69.6      71.4
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                   estimate    SE  df z.ratio p.value
##  open - shrub.ambient         -0.247 0.418 Inf -0.591  0.9764 
##  open - shrub.soil            -3.337 0.476 Inf -7.011  <.0001 
##  open - square                 1.005 0.441 Inf  2.277  0.1525 
##  open - triangle               3.195 0.517 Inf  6.179  <.0001 
##  shrub.ambient - shrub.soil   -3.090 0.551 Inf -5.612  <.0001 
##  shrub.ambient - square        1.252 0.521 Inf  2.403  0.1144 
##  shrub.ambient - triangle      3.442 0.586 Inf  5.869  <.0001 
##  shrub.soil - square           4.342 0.569 Inf  7.636  <.0001 
##  shrub.soil - triangle         6.532 0.629 Inf 10.382  <.0001 
##  square - triangle             2.190 0.603 Inf  3.629  0.0026 
## 
## P value adjustment: tukey method for comparing a family of 5 estimates
emmeans (lm.intensity, pairwise~microsite)
## $emmeans
##  microsite     emmean      SE  df asymp.LCL asymp.UCL
##  open           8.111 0.01800 Inf     8.076     8.146
##  shrub.ambient  7.395 0.04146 Inf     7.314     7.476
##  shrub.soil     3.522 0.41420 Inf     2.711     4.334
##  square         7.424 0.04371 Inf     7.338     7.510
##  triangle       7.529 0.05124 Inf     7.429     7.630
## 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                   estimate     SE  df z.ratio p.value
##  open - shrub.ambient         0.7161 0.0452 Inf 15.842  <.0001 
##  open - shrub.soil            4.5889 0.4146 Inf 11.068  <.0001 
##  open - square                0.6873 0.0473 Inf 14.539  <.0001 
##  open - triangle              0.5821 0.0543 Inf 10.718  <.0001 
##  shrub.ambient - shrub.soil   3.8728 0.4163 Inf  9.304  <.0001 
##  shrub.ambient - square      -0.0287 0.0602 Inf -0.477  0.9895 
##  shrub.ambient - triangle    -0.1340 0.0659 Inf -2.032  0.2505 
##  shrub.soil - square         -3.9015 0.4165 Inf -9.367  <.0001 
##  shrub.soil - triangle       -4.0068 0.4174 Inf -9.600  <.0001 
##  square - triangle           -0.1052 0.0674 Inf -1.562  0.5218 
## 
## Results are given on the log (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 5 estimates
#let's look at shape and blockage more closely
lm.temp.cover <- glm(temp~as.factor(microsite)* as.factor(cover.type), data = shelter.shrub.open, family="gaussian")
summary(lm.temp.cover) 
## 
## Call:
## glm(formula = temp ~ as.factor(microsite) * as.factor(cover.type), 
##     family = "gaussian", data = shelter.shrub.open)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -40.064  -16.940   -5.098   15.257   84.418  
## 
## Coefficients: (4 not defined because of singularities)
##                                                      Estimate Std. Error
## (Intercept)                                           73.6902     0.2168
## as.factor(microsite)square                            30.1200    21.5514
## as.factor(microsite)triangle                          26.5928    21.5320
## as.factor(cover.type)15                              -26.4645    21.5505
## as.factor(cover.type)50                              -30.1952    21.5505
## as.factor(cover.type)90                              -31.2449    21.5405
## as.factor(microsite)square:as.factor(cover.type)15    -2.4960     1.4509
## as.factor(microsite)triangle:as.factor(cover.type)15       NA         NA
## as.factor(microsite)square:as.factor(cover.type)50    -2.9435     1.4487
## as.factor(microsite)triangle:as.factor(cover.type)50       NA         NA
## as.factor(microsite)square:as.factor(cover.type)90         NA         NA
## as.factor(microsite)triangle:as.factor(cover.type)90       NA         NA
##                                                      t value Pr(>|t|)    
## (Intercept)                                          339.849   <2e-16 ***
## as.factor(microsite)square                             1.398   0.1623    
## as.factor(microsite)triangle                           1.235   0.2168    
## as.factor(cover.type)15                               -1.228   0.2195    
## as.factor(cover.type)50                               -1.401   0.1612    
## as.factor(cover.type)90                               -1.451   0.1469    
## as.factor(microsite)square:as.factor(cover.type)15    -1.720   0.0854 .  
## as.factor(microsite)triangle:as.factor(cover.type)15      NA       NA    
## as.factor(microsite)square:as.factor(cover.type)50    -2.032   0.0422 *  
## as.factor(microsite)triangle:as.factor(cover.type)50      NA       NA    
## as.factor(microsite)square:as.factor(cover.type)90        NA       NA    
## as.factor(microsite)triangle:as.factor(cover.type)90      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 463.5796)
## 
##     Null deviance: 7168340  on 15388  degrees of freedom
## Residual deviance: 7130318  on 15381  degrees of freedom
##   (6571 observations deleted due to missingness)
## AIC: 138155
## 
## Number of Fisher Scoring iterations: 2
anova (lm.temp.cover, test="Chisq")#all variables are significant
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: temp
## 
## Terms added sequentially (first to last)
## 
## 
##                                            Df Deviance Resid. Df
## NULL                                                       15388
## as.factor(microsite)                        2  18912.1     15386
## as.factor(cover.type)                       3  16739.0     15383
## as.factor(microsite):as.factor(cover.type)  2   2370.6     15381
##                                            Resid. Dev  Pr(>Chi)    
## NULL                                          7168340              
## as.factor(microsite)                          7149427 1.384e-09 ***
## as.factor(cover.type)                         7132688 7.104e-08 ***
## as.factor(microsite):as.factor(cover.type)    7130318   0.07755 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(lm.temp.cover, pairwise~microsite|cover.type)# square-triangle different at 15% and 90%
## $emmeans
## cover.type =  0:
##  microsite emmean     SE  df asymp.LCL asymp.UCL
##  open        73.7  0.217 Inf      73.3      74.1
##  square    nonEst     NA  NA        NA        NA
##  triangle   100.3 21.531 Inf      58.1     142.5
## 
## cover.type = 15:
##  microsite emmean     SE  df asymp.LCL asymp.UCL
##  open      nonEst     NA  NA        NA        NA
##  square      74.8  0.649 Inf      73.6      76.1
##  triangle    73.8  0.920 Inf      72.0      75.6
## 
## cover.type = 50:
##  microsite emmean     SE  df asymp.LCL asymp.UCL
##  open      nonEst     NA  NA        NA        NA
##  square      70.7  0.644 Inf      69.4      71.9
##  triangle    70.1  0.920 Inf      68.3      71.9
## 
## cover.type = 90:
##  microsite emmean     SE  df asymp.LCL asymp.UCL
##  open      nonEst     NA  NA        NA        NA
##  square      72.6  0.650 Inf      71.3      73.8
##  triangle    69.0  0.644 Inf      67.8      70.3
## 
## Confidence level used: 0.95 
## 
## $contrasts
## cover.type =  0:
##  contrast          estimate     SE  df z.ratio p.value
##  open - square       nonEst     NA  NA     NA      NA 
##  open - triangle    -26.593 21.532 Inf -1.235  0.4324 
##  square - triangle   nonEst     NA  NA     NA      NA 
## 
## cover.type = 15:
##  contrast          estimate     SE  df z.ratio p.value
##  open - square       nonEst     NA  NA     NA      NA 
##  open - triangle     nonEst     NA  NA     NA      NA 
##  square - triangle    1.031  1.126 Inf  0.916  0.6301 
## 
## cover.type = 50:
##  contrast          estimate     SE  df z.ratio p.value
##  open - square       nonEst     NA  NA     NA      NA 
##  open - triangle     nonEst     NA  NA     NA      NA 
##  square - triangle    0.584  1.123 Inf  0.520  0.8616 
## 
## cover.type = 90:
##  contrast          estimate     SE  df z.ratio p.value
##  open - square       nonEst     NA  NA     NA      NA 
##  open - triangle     nonEst     NA  NA     NA      NA 
##  square - triangle    3.527  0.915 Inf  3.853  0.0003 
## 
## P value adjustment: tukey method for comparing a family of 3 estimates
#calculating Relative Interaction Index (RII)


#turn table into wide format
data.wide.temp <- reshape (shelter.shrub.open, timevar = "microsite", v.names = "temp", direction = "wide", idvar="date")
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar,
## varying = varying, : some constant variables (month,day,hour.
## 24,hour,lat,lon,rep,intensity,pendant.id,timeblock,soil.moisture,cover.type)
## are really varying
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=triangle: first taken
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=open: first taken
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=square: first taken
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=shrub.soil: first taken
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=shrub.ambient: first taken
data.wide.intensity<- reshape(shelter.shrub.open, timevar = "microsite", v.names = "intensity", direction= "wide", idvar = "date")
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar,
## varying = varying, : some constant variables (month,day,hour.
## 24,hour,lat,lon,rep,temp,pendant.id,timeblock,soil.moisture,cover.type) are
## really varying
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=triangle: first taken
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=open: first taken
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=square: first taken
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=shrub.soil: first taken
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=shrub.ambient: first taken
#rii for all temperature data
rii.temp.shrub<- data.wide.temp %>%
  mutate(rii_calc_shrub = (temp.shrub.ambient-temp.open)/(temp.shrub.ambient + temp.open))

rii.temp.triangle<- data.wide.temp %>%
  mutate(rii_calc_triangle = (temp.triangle-temp.open)/(temp.triangle + temp.open)) 

rii.temp.square<- data.wide.temp %>%
  mutate(rii_calc_square = (temp.square-temp.open)/(temp.square+ temp.open)) 

#rii for all intensity data
rii.intensity.shrub<- data.wide.intensity %>%
  mutate(rii_calc_shrub = (intensity.shrub.ambient-intensity.open)/(intensity.shrub.ambient+ intensity.open)) #it doesn't work as well as temperature

x <- select(rii.temp.shrub, rii_calc_shrub)
y <- select(rii.temp.triangle, rii_calc_triangle)
z <- select(rii.temp.square, rii_calc_square)
rii.final.temp<-cbind(x, y, z)

rii.csv<- read.csv("C:/Users/Nargol Ghazian/Desktop/Animal-Behaviour-and-Climate-project/CH3/RII.microsite.csv")

ggplot(rii.csv, aes(Microsite, rii)) +
  geom_boxplot()+ ylab("Relative Interaction Index (RII)")+ geom_hline(yintercept=0, linetype="dashed", color = "red")+ theme_classic()

ggplot(rii.csv, aes(rii, fill = Microsite)) +
  geom_histogram() +
  geom_vline(xintercept = 0, col = 2, lty = 2) +
  labs(x = "Relative Interaction Index (RII)", y = "Frequency", fill = "") +
  scale_fill_brewer(palette = "Paired")+ theme_classic()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggline (rii.csv, x="Microsite", y= "rii", add = "mean_se", ylab = "Relative Interaction Index (RII)", xlab = "Microsite")

#GLM for RII
summary (rii.csv)
##          Microsite       rii            
##  shrub.ambient:14   Min.   :-0.0409504  
##  square       :24   1st Qu.:-0.0088791  
##  triangle     :24   Median : 0.0006209  
##                     Mean   : 0.0025064  
##                     3rd Qu.: 0.0110682  
##                     Max.   : 0.0709935
lm.rii <- glm(rii~as.factor(Microsite), data = rii.csv, family="gaussian")
summary(lm.rii) 
## 
## Call:
## glm(formula = rii ~ as.factor(Microsite), family = "gaussian", 
##     data = rii.csv)
## 
## Deviance Residuals: 
##       Min         1Q     Median         3Q        Max  
## -0.044037  -0.009890  -0.002490   0.007168   0.067907  
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)
## (Intercept)                   0.003086   0.005347   0.577    0.566
## as.factor(Microsite)square   -0.006171   0.006728  -0.917    0.363
## as.factor(Microsite)triangle  0.004673   0.006728   0.695    0.490
## 
## (Dispersion parameter for gaussian family taken to be 0.0004002093)
## 
##     Null deviance: 0.025029  on 61  degrees of freedom
## Residual deviance: 0.023612  on 59  degrees of freedom
## AIC: -304.19
## 
## Number of Fisher Scoring iterations: 2
anova (lm.rii, test="Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: rii
## 
## Terms added sequentially (first to last)
## 
## 
##                      Df  Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL                                     61   0.025029         
## as.factor(Microsite)  2 0.0014171        59   0.023612   0.1703
emmeans(lm.rii, pairwise~Microsite)#no significant difference
## $emmeans
##  Microsite       emmean      SE  df asymp.LCL asymp.UCL
##  shrub.ambient  0.00309 0.00535 Inf -0.007393   0.01357
##  square        -0.00308 0.00408 Inf -0.011088   0.00492
##  triangle       0.00776 0.00408 Inf -0.000245   0.01576
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                 estimate      SE  df z.ratio p.value
##  shrub.ambient - square    0.00617 0.00673 Inf  0.917  0.6294 
##  shrub.ambient - triangle -0.00467 0.00673 Inf -0.695  0.7667 
##  square - triangle        -0.01084 0.00578 Inf -1.878  0.1451 
## 
## P value adjustment: tukey method for comparing a family of 3 estimates
library(sjPlot)
tab_model(lm.rii)
  rii
Predictors Estimates CI p
(Intercept) 0.00 -0.01 – 0.01 0.566
Microsite [square] -0.01 -0.02 – 0.01 0.363
Microsite [triangle] 0.00 -0.01 – 0.02 0.490
Observations 62
R2 Nagelkerke 0.057 
ggplot (shelter.shrub.open, aes(as.factor(day), temp)) + geom_boxplot() + xlab("Day") + ylab ("Temperature (°F)")+ theme_classic()+ theme(axis.text=element_text(size=12))+ theme(axis.text=element_text(size=12))+facet_grid("microsite")+ stat_summary(fun.y=mean, colour="red", geom="point", shape=18, size=1,show_guide = FALSE)
## Warning: `show_guide` has been deprecated. Please use `show.legend`
## instead.

ggplot (shelter.shrub.open, aes(as.factor(day), intensity)) + geom_boxplot() + xlab("Day") + ylab ("Solar Radiation (lum/ft²)")+ theme_classic()+ theme(axis.text=element_text(size=12))+ theme(axis.text=element_text(size=12))+facet_grid("microsite")+ stat_summary(fun.y=mean, colour="red", geom="point", shape=18, size=1,show_guide = FALSE)
## Warning: `show_guide` has been deprecated. Please use `show.legend`
## instead.
## Warning: Removed 9330 rows containing non-finite values (stat_boxplot).
## Warning: Removed 9330 rows containing non-finite values (stat_summary).

#let's create a heatmap
ggplot(shelter.shrub.open, aes(as.factor(day), as.factor(microsite), fill=temp))+ geom_tile()+ xlab("Day") + ylab ("Microsite")+theme(axis.text.x = element_text(angle = 90))+ scale_fill_distiller (palette='Blues')+ theme_classic()+ labs(fill = "Temperature (°F)")

#geom_jitter plot for max and min of intensity for each microsite
ggplot(shelter.shrub.open, aes(x = as.factor(microsite), y = intensity)) +
  scale_y_log10() +
  geom_jitter(position = position_jitter(width = 0.1, height = 0), alpha = 1/4) +
  stat_summary(fun.y = min, colour = "blue", geom = "point", size = 5) +
  stat_summary(fun.y = max, colour = "red", geom = "point", size = 5)+ xlab("")+ ylab("Solar Radiation (lum/ft²)")+theme_classic()
## Warning: Removed 9330 rows containing non-finite values (stat_summary).

## Warning: Removed 9330 rows containing non-finite values (stat_summary).
## Warning: Removed 9330 rows containing missing values (geom_point).

#faceted by microsite
ggplot(shelter.shrub.open, aes(x = as.factor(day), y = intensity)) +
  scale_y_log10() +
  geom_jitter(position = position_jitter(width = 0.1, height = 0), alpha = 1/4) +
  stat_summary(fun.y = min, colour = "blue", geom = "point", size = 5) +
  stat_summary(fun.y = max, colour = "red", geom = "point", size = 5)+ xlab("Day")+ ylab("Light Intensity (lum/ft²)")+theme_classic()+ facet_grid(~microsite)+ coord_flip()+ theme(axis.text.x = element_text(angle = 90))
## Warning: Removed 9330 rows containing non-finite values (stat_summary).
## Warning: Removed 9330 rows containing non-finite values (stat_summary).
## Warning: Removed 9330 rows containing missing values (geom_point).

#MACRO-MICRO CLIMATE CONTRAST
##test between macro-climate and micro-climate

library (ggplot2)
library(sjPlot)
ggplot(macro.micro.contrast, aes((microsite), temp, fill=microsite)) + geom_boxplot() + xlab("Site") + ylab ("Temperature (°F)")+ theme_classic()+ theme(axis.text=element_text(size=12))+theme(axis.text.x = element_text(angle = 90))+ theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())+ labs(fill = "Site")+ stat_summary(fun.y=mean, colour="black", geom="point", shape=18, size=3,show_guide = FALSE)
## Warning: `show_guide` has been deprecated. Please use `show.legend`
## instead.

#weather station data are significantly differnt from shrub, open, and square.

ggline(macro.micro.contrast, x = "day", y = "temp", color = "microsite",
 add = "mean_ci", shape = "microsite", xlab = "Day", ylab = "Temperature (°F)", legend.title= "Site", legend="right")

ggline(macro.micro.contrast, x = "day", y = "intensity", color = "microsite",
 add = "mean_ci", shape = "microsite", xlab = "Day", ylab = "Solar Rdation (lum/ft²)", legend.title= "Site", legend="right")
## Warning: Removed 9330 rows containing non-finite values (stat_summary).

lm.site <- glm(temp~as.factor(microsite), data = macro.micro.contrast, family="gaussian")
summary (lm.site)
## 
## Call:
## glm(formula = temp ~ as.factor(microsite), family = "gaussian", 
##     data = macro.micro.contrast)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -39.543  -16.315   -5.344   13.475  110.218  
## 
## Coefficients:
##                                     Estimate Std. Error t value Pr(>|t|)
## (Intercept)                          73.6902     0.2193 335.970  < 2e-16
## as.factor(microsite)shrub.ambient     0.2471     0.4140   0.597   0.5506
## as.factor(microsite)shrub.soil        3.3369     0.4715   7.077 1.51e-12
## as.factor(microsite)square           -1.0049     0.4373  -2.298   0.0216
## as.factor(microsite)triangle         -3.1949     0.5122  -6.238 4.52e-10
## as.factor(microsite)weather.station  -5.4202     0.8990  -6.029 1.68e-09
##                                        
## (Intercept)                         ***
## as.factor(microsite)shrub.ambient      
## as.factor(microsite)shrub.soil      ***
## as.factor(microsite)square          *  
## as.factor(microsite)triangle        ***
## as.factor(microsite)weather.station ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 474.3471)
## 
##     Null deviance: 10784021  on 22583  degrees of freedom
## Residual deviance: 10709809  on 22578  degrees of freedom
## AIC: 203260
## 
## Number of Fisher Scoring iterations: 2
tab_model(lm.site)
  temp
Predictors Estimates CI p
(Intercept) 73.69 73.26 – 74.12 <0.001
microsite [shrub.ambient] 0.25 -0.56 – 1.06 0.551
microsite [shrub.soil] 3.34 2.41 – 4.26 <0.001
microsite [square] -1.00 -1.86 – -0.15 0.022
microsite [triangle] -3.19 -4.20 – -2.19 <0.001
microsite
[weather.station] 		-5.42 -7.18 – -3.66 <0.001
Observations 22584
R2 Nagelkerke 0.963 
library (emmeans)
anova(lm.site, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: temp
## 
## Terms added sequentially (first to last)
## 
## 
##                      Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                                 22583   10784021              
## as.factor(microsite)  5    74212     22578   10709809 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(lm.site, pairwise~microsite)
## $emmeans
##  microsite       emmean    SE  df asymp.LCL asymp.UCL
##  open              73.7 0.219 Inf      73.3      74.1
##  shrub.ambient     73.9 0.351 Inf      73.2      74.6
##  shrub.soil        77.0 0.417 Inf      76.2      77.8
##  square            72.7 0.378 Inf      71.9      73.4
##  triangle          70.5 0.463 Inf      69.6      71.4
##  weather.station   68.3 0.872 Inf      66.6      70.0
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                        estimate    SE  df z.ratio p.value
##  open - shrub.ambient              -0.247 0.414 Inf -0.597  0.9913 
##  open - shrub.soil                 -3.337 0.471 Inf -7.077  <.0001 
##  open - square                      1.005 0.437 Inf  2.298  0.1947 
##  open - triangle                    3.195 0.512 Inf  6.238  <.0001 
##  open - weather.station             5.420 0.899 Inf  6.029  <.0001 
##  shrub.ambient - shrub.soil        -3.090 0.545 Inf -5.665  <.0001 
##  shrub.ambient - square             1.252 0.516 Inf  2.426  0.1473 
##  shrub.ambient - triangle           3.442 0.581 Inf  5.925  <.0001 
##  shrub.ambient - weather.station    5.667 0.940 Inf  6.030  <.0001 
##  shrub.soil - square                4.342 0.563 Inf  7.708  <.0001 
##  shrub.soil - triangle              6.532 0.623 Inf 10.480  <.0001 
##  shrub.soil - weather.station       8.757 0.967 Inf  9.059  <.0001 
##  square - triangle                  2.190 0.598 Inf  3.664  0.0034 
##  square - weather.station           4.415 0.950 Inf  4.646  <.0001 
##  triangle - weather.station         2.225 0.987 Inf  2.254  0.2131 
## 
## P value adjustment: tukey method for comparing a family of 6 estimates
lm.site.sun <- glm(intensity~as.factor(microsite), data = macro.micro.contrast, family="gaussian")
summary (lm.site.sun)
## 
## Call:
## glm(formula = intensity ~ as.factor(microsite), family = "gaussian", 
##     data = macro.micro.contrast)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -16478   -2436   -1303     124   47091  
## 
## Coefficients:
##                                     Estimate Std. Error t value Pr(>|t|)
## (Intercept)                          3331.50      80.51  41.382  < 2e-16
## as.factor(microsite)shrub.ambient   -1703.47     152.58 -11.165  < 2e-16
## as.factor(microsite)shrub.soil      -3297.64     203.37 -16.215  < 2e-16
## as.factor(microsite)square          -1656.02     160.31 -10.330  < 2e-16
## as.factor(microsite)triangle        -1470.12     189.26  -7.768 8.58e-15
## as.factor(microsite)weather.station 13146.49     260.91  50.387  < 2e-16
##                                        
## (Intercept)                         ***
## as.factor(microsite)shrub.ambient   ***
## as.factor(microsite)shrub.soil      ***
## as.factor(microsite)square          ***
## as.factor(microsite)triangle        ***
## as.factor(microsite)weather.station ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 38433720)
## 
##     Null deviance: 6.4238e+11  on 13253  degrees of freedom
## Residual deviance: 5.0917e+11  on 13248  degrees of freedom
##   (9330 observations deleted due to missingness)
## AIC: 269095
## 
## Number of Fisher Scoring iterations: 2
tab_model(lm.site.sun)
  intensity
Predictors Estimates CI p
(Intercept) 3331.50 3173.71 – 3489.29 <0.001
microsite [shrub.ambient] -1703.47 -2002.52 – -1404.43 <0.001
microsite [shrub.soil] -3297.64 -3696.23 – -2899.05 <0.001
microsite [square] -1656.02 -1970.22 – -1341.83 <0.001
microsite [triangle] -1470.12 -1841.06 – -1099.17 <0.001
microsite
[weather.station] 		13146.49 12635.12 – 13657.86 <0.001
Observations 13254
R2 Nagelkerke 1.000 
library (emmeans)
anova(lm.site.sun, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: intensity
## 
## Terms added sequentially (first to last)
## 
## 
##                      Df   Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                                   13253 6.4238e+11              
## as.factor(microsite)  5 1.3321e+11     13248 5.0917e+11 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(lm.site.sun, pairwise~microsite)
## $emmeans
##  microsite        emmean    SE  df asymp.LCL asymp.UCL
##  open             3331.5  80.5 Inf      3174      3489
##  shrub.ambient    1628.0 129.6 Inf      1374      1882
##  shrub.soil         33.9 186.8 Inf      -332       400
##  square           1675.5 138.6 Inf      1404      1947
##  triangle         1861.4 171.3 Inf      1526      2197
##  weather.station 16478.0 248.2 Inf     15992     16964
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                        estimate  SE  df z.ratio p.value
##  open - shrub.ambient              1703.5 153 Inf  11.165 <.0001 
##  open - shrub.soil                 3297.6 203 Inf  16.215 <.0001 
##  open - square                     1656.0 160 Inf  10.330 <.0001 
##  open - triangle                   1470.1 189 Inf   7.768 <.0001 
##  open - weather.station          -13146.5 261 Inf -50.387 <.0001 
##  shrub.ambient - shrub.soil        1594.2 227 Inf   7.013 <.0001 
##  shrub.ambient - square             -47.5 190 Inf  -0.250 0.9999 
##  shrub.ambient - triangle          -233.4 215 Inf  -1.086 0.8870 
##  shrub.ambient - weather.station -14850.0 280 Inf -53.039 <.0001 
##  shrub.soil - square              -1641.6 233 Inf  -7.058 <.0001 
##  shrub.soil - triangle            -1827.5 253 Inf  -7.212 <.0001 
##  shrub.soil - weather.station    -16444.1 311 Inf -52.944 <.0001 
##  square - triangle                 -185.9 220 Inf  -0.844 0.9593 
##  square - weather.station        -14802.5 284 Inf -52.072 <.0001 
##  triangle - weather.station      -14616.6 302 Inf -48.472 <.0001 
## 
## P value adjustment: tukey method for comparing a family of 6 estimates
ggplot(macro.micro.contrast, aes(temp, fill = microsite)) +
         geom_histogram(binwidth = 5) +
  scale_fill_brewer(palette = "Set1")+   labs(fill = "", x = "Temperature (°F)", y = "Frequency")+theme_classic()+ theme(axis.text=element_text(size=12))+ stat_central_tendency(aes(color = microsite), type = "mean", color="green", linetype = 2)+ stat_central_tendency(type = "median", color = "blue", linetype = 2)+ facet_grid(~microsite)+theme(legend.position = "none")

ggplot(macro.micro.contrast, aes((day), temp, color=microsite)) + geom_smooth()+ xlab("Day") + ylab ("Temperature (°F)")+ theme_classic()+ theme(axis.text=element_text(size=12))+stat_summary(fun.y=max, geom="point", size=2, aes(shape = microsite))+ facet_grid("microsite")+ theme(legend.position = "none")
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

ggplot(macro.micro.contrast, aes((day), intensity, color=microsite)) + geom_smooth()+ xlab("Day") + ylab ("Solar Rdation (lum/ft²)")+ theme_classic()+ theme(axis.text=element_text(size=12))+stat_summary(fun.y=max, geom="point", size=2, aes(shape = microsite))+ facet_grid("microsite")+ theme(legend.position = "none")
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## Warning: Removed 9330 rows containing non-finite values (stat_smooth).
## Warning: Removed 9330 rows containing non-finite values (stat_summary).

library(car)
## Warning: package 'car' was built under R version 3.6.3
## Loading required package: carData
## Warning: package 'carData' was built under R version 3.6.3
## Registered S3 methods overwritten by 'car':
##   method                          from
##   influence.merMod                lme4
##   cooks.distance.influence.merMod lme4
##   dfbeta.influence.merMod         lme4
##   dfbetas.influence.merMod        lme4
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
## The following object is masked from 'package:purrr':
## 
##     some
leveneTest(temp ~ microsite, macro.micro.contrast)
## Levene's Test for Homogeneity of Variance (center = median)
##          Df F value    Pr(>F)    
## group     5  60.096 < 2.2e-16 ***
##       22578                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
leveneTest(intensity ~microsite, macro.micro.contrast)
## Levene's Test for Homogeneity of Variance (center = median)
##          Df F value    Pr(>F)    
## group     5  815.31 < 2.2e-16 ***
##       13248                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#variation in sunlight and intensity are significant


wide.clim <- reshape(macro.micro.contrast, v.names="temp", timevar="microsite", idvar=c("temp", "microsite"),
        direction="wide")
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar,
## varying = varying, : some constant variables (date,month,day,hour.
## 24,hour,lat,lon,rep,intensity,pendant.id,timeblock,soil.moisture,cover.type)
## are really varying
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=triangle: first taken
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=open: first taken
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=square: first taken
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=shrub.soil: first taken
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=shrub.ambient: first taken
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=weather.station: first taken
summary(wide.clim)
##       site               sensor           date           year     
##  panoche:2671   hobo.pendant:2429   05/27/19: 168   Min.   :2019  
##                 satellite   : 242   05/24/19: 133   1st Qu.:2019  
##                                     05/21/19: 130   Median :2019  
##                                     05/20/19: 126   Mean   :2019  
##                                     05/25/19: 126   3rd Qu.:2019  
##                                     05/29/19: 123   Max.   :2019  
##                                     (Other) :1865                 
##      month            day           hour.24              hour     
##  Min.   :5.000   Min.   : 1.0   10:46:02:  22              : 242  
##  1st Qu.:5.000   1st Qu.: 7.0   8:22:41 :  22   10:46:02 AM:  22  
##  Median :5.000   Median :21.0   1:38:16 :  21   8:22:41 AM :  22  
##  Mean   :5.364   Mean   :17.6   10:22:04:  21   1:38:16 AM :  21  
##  3rd Qu.:6.000   3rd Qu.:26.0   10:38:01:  20   10:22:04 AM:  21  
##  Max.   :6.000   Max.   :31.0   12:38:43:  20   1:24:11 PM :  20  
##                                 (Other) :2545   (Other)    :2323  
##       lat             lon              rep           intensity    
##  Min.   :36.39   Min.   :-120.9   Min.   : 1.000   Min.   :    0  
##  1st Qu.:36.69   1st Qu.:-120.8   1st Qu.: 1.000   1st Qu.:   96  
##  Median :36.69   Median :-120.8   Median : 2.000   Median :  768  
##  Mean   :36.69   Mean   :-120.8   Mean   : 2.404   Mean   : 5697  
##  3rd Qu.:36.70   3rd Qu.:-120.8   3rd Qu.: 3.000   3rd Qu.: 4864  
##  Max.   :36.70   Max.   :-120.8   Max.   :12.000   Max.   :63505  
##  NA's   :242     NA's   :242      NA's   :242      NA's   :727    
##    pendant.id           timeblock    soil.moisture      cover.type   
##  Min.   : 1156629            :1024   Min.   :  1.00   Min.   : 0.00  
##  1st Qu.: 1174238   morning  : 714   1st Qu.: 23.00   1st Qu.: 0.00  
##  Median :20567072   evening  : 563   Median : 29.40   Median :15.00  
##  Mean   :13343324   afternoon: 128   Mean   : 82.71   Mean   :21.51  
##  3rd Qu.:20567073   0        :  50   3rd Qu.: 31.80   3rd Qu.:15.00  
##  Max.   :20619223   0.0001   :   3   Max.   :998.00   Max.   :90.00  
##  NA's   :242        (Other)  : 189   NA's   :1074     NA's   :1266   
##  temp.triangle      temp.open       temp.square     temp.shrub.soil 
##  Min.   : 36.53   Min.   : 34.20   Min.   : 34.79   Min.   : 45.10  
##  1st Qu.: 56.84   1st Qu.: 58.22   1st Qu.: 58.56   1st Qu.: 66.62  
##  Median : 74.62   Median : 80.21   Median : 79.15   Median : 87.55  
##  Mean   : 75.37   Mean   : 81.84   Mean   : 80.63   Mean   : 91.82  
##  3rd Qu.: 93.21   3rd Qu.:104.23   3rd Qu.:101.38   3rd Qu.:113.47  
##  Max.   :124.08   Max.   :156.98   Max.   :156.98   Max.   :187.25  
##  NA's   :2256     NA's   :2160     NA's   :2192     NA's   :2192    
##  temp.shrub.ambient temp.weather.station
##  Min.   : 34.39     Min.   : 43.88      
##  1st Qu.: 60.54     1st Qu.: 58.87      
##  Median : 84.11     Median : 69.89      
##  Mean   : 87.14     Mean   : 70.90      
##  3rd Qu.:110.92     3rd Qu.: 82.72      
##  Max.   :162.72     Max.   :100.58      
##  NA's   :2126       NA's   :2429
wide.clim.sun <- reshape(macro.micro.contrast, v.names="intensity", timevar="microsite", idvar=c("temp", "microsite"), direction="wide")
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar,
## varying = varying, : some constant variables (date,month,day,hour.
## 24,hour,lat,lon,rep,pendant.id,timeblock,soil.moisture,cover.type) are
## really varying
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=triangle: first taken
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=open: first taken
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=square: first taken
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=shrub.soil: first taken
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=shrub.ambient: first taken
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=weather.station: first taken
summary (wide.clim.sun)
##       site               sensor           date           year     
##  panoche:2671   hobo.pendant:2429   05/27/19: 168   Min.   :2019  
##                 satellite   : 242   05/24/19: 133   1st Qu.:2019  
##                                     05/21/19: 130   Median :2019  
##                                     05/20/19: 126   Mean   :2019  
##                                     05/25/19: 126   3rd Qu.:2019  
##                                     05/29/19: 123   Max.   :2019  
##                                     (Other) :1865                 
##      month            day           hour.24              hour     
##  Min.   :5.000   Min.   : 1.0   10:46:02:  22              : 242  
##  1st Qu.:5.000   1st Qu.: 7.0   8:22:41 :  22   10:46:02 AM:  22  
##  Median :5.000   Median :21.0   1:38:16 :  21   8:22:41 AM :  22  
##  Mean   :5.364   Mean   :17.6   10:22:04:  21   1:38:16 AM :  21  
##  3rd Qu.:6.000   3rd Qu.:26.0   10:38:01:  20   10:22:04 AM:  21  
##  Max.   :6.000   Max.   :31.0   12:38:43:  20   1:24:11 PM :  20  
##                                 (Other) :2545   (Other)    :2323  
##       lat             lon              rep              temp       
##  Min.   :36.39   Min.   :-120.9   Min.   : 1.000   Min.   : 34.20  
##  1st Qu.:36.69   1st Qu.:-120.8   1st Qu.: 1.000   1st Qu.: 59.88  
##  Median :36.69   Median :-120.8   Median : 2.000   Median : 79.50  
##  Mean   :36.69   Mean   :-120.8   Mean   : 2.404   Mean   : 82.50  
##  3rd Qu.:36.70   3rd Qu.:-120.8   3rd Qu.: 3.000   3rd Qu.:101.69  
##  Max.   :36.70   Max.   :-120.8   Max.   :12.000   Max.   :187.25  
##  NA's   :242     NA's   :242      NA's   :242                      
##    pendant.id           timeblock    soil.moisture      cover.type   
##  Min.   : 1156629            :1024   Min.   :  1.00   Min.   : 0.00  
##  1st Qu.: 1174238   morning  : 714   1st Qu.: 23.00   1st Qu.: 0.00  
##  Median :20567072   evening  : 563   Median : 29.40   Median :15.00  
##  Mean   :13343324   afternoon: 128   Mean   : 82.71   Mean   :21.51  
##  3rd Qu.:20567073   0        :  50   3rd Qu.: 31.80   3rd Qu.:15.00  
##  Max.   :20619223   0.0001   :   3   Max.   :998.00   Max.   :90.00  
##  NA's   :242        (Other)  : 189   NA's   :1074     NA's   :1266   
##  intensity.triangle intensity.open  intensity.square intensity.shrub.soil
##  Min.   :    2      Min.   :    1   Min.   :    1    Min.   :   1.00     
##  1st Qu.:  384      1st Qu.:  468   1st Qu.:  375    1st Qu.:   4.50     
##  Median : 1184      Median : 1752   Median : 1184    Median :  16.00     
##  Mean   : 3938      Mean   : 5631   Mean   : 2376    Mean   :  53.68     
##  3rd Qu.: 5568      3rd Qu.: 9216   3rd Qu.: 3712    3rd Qu.:  29.00     
##  Max.   :19456      Max.   :27648   Max.   :24576    Max.   :1152.00     
##  NA's   :2381       NA's   :2299    NA's   :2329     NA's   :2380        
##  intensity.shrub.ambient intensity.weather.station
##  Min.   :    1           Min.   :    0.0          
##  1st Qu.:  248           1st Qu.:  334.1          
##  Median :  864           Median :17467.0          
##  Mean   : 2476           Mean   :24799.5          
##  3rd Qu.: 3968           3rd Qu.:48678.7          
##  Max.   :19456           Max.   :63504.9          
##  NA's   :2264            NA's   :2429
#Create microsite map for appendix
library(ggmap)
## Warning: package 'ggmap' was built under R version 3.6.3
## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.
## Please cite ggmap if you use it! See citation("ggmap") for details.
## 
## Attaching package: 'ggmap'
## The following object is masked from 'package:magrittr':
## 
##     inset
register_google(key="AIzaSyBpfKtYrkYVS3LEJSjV1cIHeYrxJPsPX4U")
panoche1 <- get_map(location = c(lon = -120.7932, lat = 36.69363), zoom = 6, maptype = "satellite")
## Source : https://maps.googleapis.com/maps/api/staticmap?center=36.69363,-120.7932&zoom=6&size=640x640&scale=2&maptype=satellite&language=en-EN&key=xxx
panoche2 <- get_map(location = c(lon = -120.7932, lat = 36.69363), zoom = 12, maptype = "terrain")
## Source : https://maps.googleapis.com/maps/api/staticmap?center=36.69363,-120.7932&zoom=12&size=640x640&scale=2&maptype=terrain&language=en-EN&key=xxx
lat_long <- read.csv("C:/Users/Nargol Ghazian/Desktop/Animal-Behaviour-and-Climate-project/CH3/Shrub_Shelter_Open lat long.csv")
na.omit(lat_long)
##    Shelter.ID      Lat      Long Microsite
## 1           1 36.69363 -120.7932  Triangle
## 2           2 36.69364 -120.7933    Square
## 3           3 36.69355 -120.7931    Square
## 4           4 36.69349 -120.7932  Triangle
## 5           5 36.69349 -120.7931  Triangle
## 6           6 36.39342 -120.7931    Square
## 7           7 36.69394 -120.7930    Square
## 8           8 36.69397 -120.7929  Triangle
## 9           9 36.69401 -120.7928    Square
## 10         10 36.69400 -120.7930  Triangle
## 11         11 36.69405 -120.7930    Square
## 12         12 36.69408 -120.7930  Triangle
## 13          1 36.69534 -120.7970     Shrub
## 14          2 36.69532 -120.7972     Shrub
## 15          3 36.69534 -120.7972     Shrub
## 16          4 36.69606 -120.7967     Shrub
## 17          5 36.69593 -120.7967     Shrub
## 18          6 36.69596 -120.7968     Shrub
## 19          7 36.69589 -120.7980     Shrub
## 20          1 36.69529 -120.7970      Open
## 21          2 36.69533 -120.7971      Open
## 22          3 36.69529 -120.7923      Open
## 23          4 36.69596 -120.7966      Open
## 24          5 36.69587 -120.7967      Open
## 25          6 36.69601 -120.7967      Open
## 26          7 36.69590 -120.7981      Open
ggmap(panoche2) +
  geom_point(data= lat_long, aes(x=Long, y=Lat, color = Microsite), alpha = 6/10, size =2, show.legend = TRUE) +
  labs(title = "Map of Shelters and Shrubs", x = "Longitude", y = "Latitude")
## Warning: Removed 1 rows containing missing values (geom_point).

ggmap(panoche1) +
  geom_point(data= lat_long, aes(x=Long, y=Lat, color = Microsite), alpha = 6/10, size =2, show.legend = TRUE) +
  labs(title = "Map of Califronia", x = "Longitude", y = "Latitude")

#Very general arial view of microsites
#CONCLUSIONS

##1.Temperature and sunlight intensity are positively related.
##2.Shelters function similar to shrubs, but square emulates shrub slightly better. 
##3.Triangle at 90% blockage is best at lowering temperature. 
##4.The open experienced the most variation. 
##5. Temperature from Weather station data were significantly lower than square, shrub, and the open- weather-station data is nor ecologically-relevant.